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ABSTRACT 

Using a suite of self-gravitating, coUisionless A^-body models, we systematically explore a parameter 
space relevant to the onset and behavior of the radial orbit instability (ROI), whose strength is 
measured by the systemic axis ratios of the models. We show that a combination of two initial 
conditions, namely the velocity anisotropy and the virial ratio, determines whether a system will 
undergo ROI and exactly how triaxial the system will become. A third initial condition, the radial 
shape of the density profile, plays a smaller, but noticeable role. Regarding the dynamical development 
of the ROI, the instability a) begins after systems collapse to their most compact configuration and 
b) evolves fastest when a majority of the particles have radially anisotropic orbits while there is a lack 
of centrally-concentrated isotropic orbits. We argue that this is further evidence that self-reinforcing 
torques are the key to the onset of the ROI. Our findings support the idea that a separate orbit 
instability plays a role in halting the ROI. 

Subject headings: galaxics:structure — galaxics:kincniatics and dynamics 



1. INTRODUCTION 

A^-body simulations of self-gravitating, coUisionless 
systems of massive particles have become a standard 
tool used to investigate the behavior of galactic systems. 
Such particle aggregations can represent the visible, stel- 
lar components of galaxies as well as the dark matter 
that is hypothesized to surround individual galaxies (ha- 
los) . Since the dark matter halos of galaxies are thought 
to contain the majority of the mass in a galactic system, 
the simplest approximation to modeling a real galaxy is 
to ignore the stars and focus on the dark matter. 

This approach has led to cosmological simulations that 
produce cosmic webs of dark matter, where hierarchi- 
cal merging creates sub-galactic, galactic, and cluster 
sized systems. Mo re sophisticated mod els (like the Mil- 
lenium Simulation. ISpringel et aLll2005[ ) also include gas 
and stars, but the dark matter remains the dominant 
gravitational component. Simulations like these have 
found a number of interesting properties of dark mat- 
ter systems. The density profiles of dark matter ha- 
los t end to have shapes described by a simple relation 
(e.Q.. lNavarro "alll2004[ ) . Generically, density increases 
with decreasing radius, but the increase is larger near 
the edges of systems. The overall shape of the den- 
sity profiles appears to be a robust outcome o f dark 
matter halo e v olution (iLe Delliou fc Henri kserJ l2003l : 
iLu et al\ 120061: iGonzalez-Casado et a l. 200^^^ How- 
ever, the inclusion of gaseous baryonic material can 
impact the innermost regions of dark matter halos 
(|Kazantzidis. Zentner. fc Kravtsovll2006[ ). possibly play- 
ing an important role in halo evolution. 

Our overall goal is to better understand the physics in- 
volved in the evolution of dark matter halos. Such under- 
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standing should explain, for example, the shapes of halo 
density profiles and why the coarse-g rained phase-space 
densi t y is single power-law in radius (jTavlor fc Navarrol 
I2OOII: iBarnes et all |2001D . In this work, we focus 
on models that involve, to some degree, the well- 
disc ussed radial orbit instability ( ROI) (for exam- 
ple. iPolvachenko fc ShukhmanI 1198X1: Ivan Albadal | 1982|: 



Mcrritt fc Aguilai' 'igSSl IPalmcr & Papaloizou' '1987 
Katz 1991: Saha 1991; .C in_cgtta, Nunez, fc Muzzio 199e 
Theis fc SpurzemI Il999l: IHuss. Jain, fc Steinmetd Il99£ 



Barries et al. 20051: MacMillan Widrow. fc HenrikserJ 
2006HBoilv fc Athanassoulall2006HBellovarv et aLll2008l ). 

To this end, we have constructed a set of A^-body mod- 
els with a variety of defining characteristics to represent 
dark matter halos. These models have been evolved in 
time, and the evolutions have been analyzed to find com- 
mon trends in behavior. 

The ROI can be briefly summarized as follows. Initally 
spherical A^-body systems made up of particles that have 
predominantly radial motions do not remain spherical. 
In general, the overall shapes of these systems become 
prolate spheroidal or triaxial (but on time scales of or- 
der the relaxa tion time, they can bec ome more spheri- 
cal again, see iTheis fc SpurzemI Il999t) . The instability 
appears in models with and without Hubble expansion. 
As this instability can drive changes in the density and 
velocity distributions in dark matter systems, we are in- 
terested in furthering the description and understanding 
of its mechanisms. 

We will present relationships between initial conditions 
of A^-body models and the behavior of the ROI. In par- 
ticular, we show that neither global measures of velocity 
anisotropy nor virial ratio alone can be used to predict 
the onset of the ROI. Instead, we find that combining 
these measures provides a determination of the end re- 
sult of the ROI, specifically the overall shape of the sys- 
tem. We also argue that the results of our simulations 
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su pport the idea that another orbit instability, detailed 
by I Adams et al\ ()2007l ). plays an important role in halt- 
ing the ROI. However, we point out that our simulations 
cannot speak to any possible role gas could play in stop- 
ping the ROI. 

As a first step in gathering information about the 
physics of collisionless, self-gravitating systems, we inves- 
tigate models that do not include Hubble expansion. As 
we will discuss, these non-cosmological systems contain 
complexities that would be compounded by the inclu- 
sion of initial radial expansion. Section [2] begins with an 
explanation of our model creation and evolution meth- 
ods. Our analysis routines and testing criteria are also 
described in this section. The crux of this work begins 
with § 131 a discussion of the global behavior of the ROI in 
our models. Section 2] deals with our observations of the 
onset and early development of the ROI. Executive sum- 
maries of these topics can be found in Sections l3.3[ 14.1.31 
and 14.2.41 A synopsis of this work and our conclusions 
are presented in § l5l 

2. METHODS & TESTING 
2.1. Initial Conditions 

We create A^-body systems with three different initial 
density profiles; p ~ po (constant), p oc r^^ (cuspy), and 
p cc (Gaussian). Particles are distributed randomly 
within a sphere of radius R = 1 according to the den- 
sity distribution so that the total mass M = 1. Both R 
and M are discussed in dimensionless code units where 
G = 1; fixing a physical value of M and R would al- 
low us to also transform any other quantities to physical 
units. The cuspy and Gaussian profiles involve scale- 
lengths; the density at the scalelength is half its central 
value for the cuspy profile and times its central value 
for the Gaussian profile. The cuspy scalelength is 10~^ 
and the Gaussian scalelength is l/\/5- The density pro- 
file of a cuspy system is shown in Figure [TJi. For most of 
our simulations, the particles have equal masses, but we 
have also created models with particles of different mass 
(see ? l2.3p . For these multi-mass models, the fractions of 
particles with each mass is determined so that the total 
mass of the system remains M = 1 . 

With the initial positions of the particles fixed, we cal- 
culate the potential energy of a system and adopt a value 
for the initial virial ratio Qo = to calculate a 

scale speed that is used to assign individual particle ve- 
locities. The virial ratio changes during the evolution 
of the systems, and we denote its value at times after 
t = simply as Q. Velocities are assigned to particles 
based on an assumed radial velocity anisotropy profile, 
where (3=1 — cr^an/^o'^ad ^^e velocity anisotropy and 
(T denotes velocity dispersion. We utili ze the anisotropy 
profile discussed in lBarnes et al\ (|2007l ). 

/3(^) = i(/3high-Aow)[l + tanh(TOlogr/ra)]+/?iow, (1) 

where /3iow is the anisotropy value at r = 0, /3high is the 
value at large r, m controls the transition between these 
values, and is the anisotropy radius, at which (3 = 0.5. 
The flexibility of this relation allows us to easily create 
completely isotropic models (/3 = 0), completely radially 
anisotropic models (/3 = 1), or systems with intermedi- 
ate levels of radial anisotropy. With fixed values of /3iow 



and /3higii, the amount of anisotropy is controlled mainly 
by the value of ra, but the value of the slope m also con- 
tributes. For the work discussed here, we have fixed the 
slope value m = 7 to provide a relatively rapid change 
between (3\ow and /3high- Note that a given anisotropy 
profile will provide different amounts of mass with radi- 
ally anisotropic velocities depending on the density pro- 
file. An example /3(r) with /3iow = 0, /3high = 1, and 
Ta = 0.5 is plotted in Figure [1)3. 

For a chosen anisotropy profile, we proceed to assign 
individual velocities to particles. A particle's speed can 
be taken from a distribution about the scale speed. In 
isotropic regions, particles have velocity directions that 
uniformly cover 47r steradians. Radial anisotropy is guar- 
anteed by limiting the velocity vector to lie within cones 
centered on a particle with axes along the radial direc- 
tion (one cone opens towards the center, the other opens 
outward). The opening angle of the cone correlates to 
the desired value of [3. For small /3, the angle is close to 
tt; when /3 « 1, the angle is nearly 0. The velocity vec- 
tor for a specific particle is chosen randomly within the 
allowed range of directions. Since the vector has definite 
limits placed on its direction, we refer to this as the 'hard- 
edge' setup. Panels c and d of Figure [H show the average 
component velocities and dispersions, respectively, for a 
model with the same anisotropy parameters as above and 
Qo = 0.5. The average component velocities for all shells 
in all models are initially zero. 

Of course, this is not the only route to initially 
anisotropic systems. We have also created models with 
initial conditions that contain combinations of isotropic 
and radially anisotropic orbits at all radii, in contrast 
to the 'hard-edge' setup. The evolutions of these 'soft- 
edge' models is very similar to those using the 'hard-edge' 
setup. All initially anisotropic models subsequently dis- 
cussed use the 'hard-edge' setup. 

2.2. Evolution & Analysis 

We have used the direct A^-body integration code 
NB0DY2 and the treecode Gadget-2 to evolve our en- 
sembles of particles. Readers interested in the details of 
NB0I )Y2 and Gadget -2 al gorithms are enc ouraged to ex- 
amine (AarsitB (|2001l ) and ISpringell l|2005D . respectively. 
Any two particles i and j interact through softened forces 
of the form, 



F,, = -G 
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where r.y is the distance between the particles and e is 
the softening length. NB0DY2 allows particle-particle 
merging, but we do not implement that option in our 
models. The time interval in these models is the crossing 
time determined by the initial density profile. A typi- 
cal run evolves the system for approximately 20 inital 
crossing times Across- The actual length of an evolution 
depends on the strength of the collapse, with stronger 
collapses running for fewer crossing times as close particle 
encounters decrease timestcp sizes. However, all simula- 
tions have been run long enough for the system to come 
to virial equilibrium, and most have reached mechanical 
equilibrium throughout. Some models have an outer re- 
gion containing w 5% of the total mass that continues to 
expand throughout the simulations. NB0DY2 includes 
a routine for removing particles that become unbound 
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from the system, however Gadget-2 does not. The direct 
nature of NB0DY2 Hmits the particle numbers that we 
can reasonably investigate (in this work iV = 10''). As a 
treecode, Gadget-2 allows us to evolve models with an or- 
der of magnitude more particles {N = 10^) in roughly the 
same amount of time. Aside from the ability to remove 
escapers with NB0DY2, the biggest difference between 
the two codes is that energy conservation is much better 
with NB0DY2 (but adequate with Gadget-2). 

Output from both codes occurs at fixed intervals deter- 
mined by the initial crossing time of the system fcross ■ We 
have modified the NB0DY2 ouput process somewhat to 
keep track of different variables than originally intended. 
We create "snapshots" of the system at the output in- 
tervals; numbers of particles, postions, velocities, ener- 
gies, numbers and strengths of coUsions, etc. Gadget-2 
output files containing particle positions and velocities 
are straightforwardly generated. Once a model has com- 
pleted its evolution, the snapshots are analyzed to deter- 
mine the evolution of quantities such as density, velocity 
dispersion and anisotropy, and axis ratios. 

Our analysis tools generally come in two groups, spher- 
ical and ellipsoidal, but we utilize only spherical tools in 
this work. Our analyses place particles in shells that 
each contain a fraction of the total system mass. For 
NB0DY2 models, each shell contains 5% of the bound 
mass, while Gadget-2 shells contain 1% of the total mass. 
These shells are used to determine average densities, av- 
erage velocites and dispersions, etc. that give us radial 
profiles of these quantities at any point during the run. 
Additionally, two sets of systemic axis ratios are deter- 
mined. These axis ratios relate the length of the longest 
axis of a system a to the lengths of the intermediate b and 
short c axes. One pair of (^,f ) values is created using 
positions of particles that make up the innermost 95% of 
the mass in a halo. The outer 5% of the mass has been 
discounted as it can often reach distances greater than 
10 half-mass radii. Particles at such large distances can 
allow small amounts of mass to severely distort the axis 
ratios from values appropriate for the more centrally con- 
centrated mass of the system. The other axis ratio pair 
gives the shape of the innermost 80% of the mass. Com- 
paring the behaviors of these two sets of axis ratios gives 
us leverage to understand which part of the halo (inner 
or outer) is dominating the overall shape. We note that 
since Gadget-2 docs not remove unbound particles, 95% 
mass axis ratios from Gadget-evolved systems tend to be 
rounder than those derived from NB0DY2 models. 

We also track velocity quantities throughout the evo- 
lutions. Average spherical component velocities (vr, v^, 
vg) and velocity dispersions ((7^, cr^, ae) arc calculated 
for the particles in a given shell. With the dispersions, 
we calculate the anisotropy parameter P for each shell 
(effectively creating a radial profile). We divide the 
anisotropy categories as follows: if the anisotropy of a 
shell (3 > 0.5, it is radially anisotropic; isotropic shells 
have -1.0 < (3 < 0.5, and (3 < -1.0 for a shell to be 
tangentially anisotropic. We have chosen these bound- 
ary values since /? = 0.5 occurs when the radial disper- 
sion is double that for cither of the tangential directions 
(fTr — 2(70 = 2(70). Similarly, f3 ~ —1.0 when each of 
the tangential dispersions are twice the radial dispersion 
(tT0 = ao ~ 2ar). 



We keep track of the fractions of mass with 
isotropic {^i), radially anisotropic (/ir), and tangentially 
anisotropic (fit) orbits by adding the particle numbers of 
shells with appropriate (3 values. We also calculate the 
fraction of centrally-concentrated mass that is isotropic 
fii^c, where a centrally-concentrated shell must be part 
of the innermost 50% of the mass of a system. The mass 
fractions are related through fi^ + f^t + l^-i = 1, and while 
Hi can be as large as 1, /i^^c has a maximum value of 0.50. 

2.3. Testing 

We have tested the robustness of our NB0DY2 results 
by varying two important parameters, softening length 
and number of particles. The NB0DY2 softening length 
has been varied by two orders of magnitude (e = 5 x 10^^ 
to e = 5 X 10~^). Pairs of models with identical ini- 
tial conditions and parameter sets, save softening length, 
show that softenings in this range lead to essentially iden- 
tical evolutions. Gadgct-2 simulations use e = 1 x 10"'*. 
We note that thes e values of e are a pproximately the 
r ange suggested bv iPower et "^(l200l . 

iBoilv. Athanassoula. k, Kroupal (|2002f ) and 
iBoilv k, Athanassoulal (|2006ft argue that a conser- 
vative estimate of iV « 10^ is required to accurately 
follow non-spherical collapses. To further test the 
impact of resolution on the quantities relevant for this 
work, we have also re-evolved 3 versions of = 10^ 
models using TV = 5 x 10'' with NB0DY2. These 
specific models were chosen to represent the range of 
initial density profiles (constant to Gaussian), velocity 
anisotropy (isotropic to radially anisotropic), and kinetic 
temperatures (very cold to warm). The A^ = 5 x 10^ 
models behave very similarly to the A^ = 10"* models. 
We also ran a single A^ = 10^ run with NB0DY2 as an 
extreme test. The most important quantities for this 
work, overall strength of the ROI in terms of axis ratios, 
rate of ROI onset, and velocity dispersion evolution, do 
not appear to change dramatically with increased N . 
Most importantly, the Gadget-2 models with A^ ~ 10^ 
do not show any coordinated or significant deviations in 
behavior from their NB0DY2 counterparts. 

The softening value plays a key role in the degree of 
collisionality of the evolutions. If e is too small, close 
encounters between particles give rise to unwanted two- 
body effects. When e is too large, the fundamental inter- 
action between particles ceases to be Newtonian grav- 
ity. We have tested the collisionality of our halos in 
two ways, using NB0DY2. First, we have tracked the 
changes in particle velocities At7 throughout our sim- 
ulations. Strong two-body collisions lead to values of 
|Au|/u w 1. We find that only a small fraction of the 
particles in our models undergo such strong collisions. 
Second, because collisional effects would drive mass seg- 
regation on the two-body relaxation timescale, we have 
evolved multi-mass models where particles are divided 
into two groups with different masses (in our tests, the 
mass ratio is 3). Two multi-mass models with soften- 
ings e = 5 X 10~^ and e = 5 x 10~^, respectively, have 
been evolved for approximately one two-body relaxation 
time (« 140 tcross)- The evolutions show virtually iden- 
tical behaviors; in each run, the half-mass radius of the 
heavier particles quickly reaches a value that remains 
virtually constant up to one relaxation time. The lack 
of mass segregation indicates that our simulations within 



4 



Barnes et al. 



this softening range are coUisionlcss for at least 100 cross- 
ing times. 

These tests demonstrate that the essential physics of 
collisionless, self-gravitating collapses is present in our 
simulations. Of course, larger particle numbers would 
allow higher resolution, but our goal is not to push the 
boundaries of numerical simulation. These straightfor- 
ward simulations allow us to focus on understanding the 
physics of the radial orbit instability by investigating a 
large variety of initial conditions. 

3. GLOBAL BEHAVIOR OF THE RADIAL ORBIT 
INSTABILITY 

The most common criteria for determining whether or 
not the radial orbit instability has occurred involves the 
overall shape of the system. If a system remains more-or- 
less spherical throughout its evolution, the ROI has not 
developed. When a system changes from an initial spher- 
ical shape to a final spheroidal or triaxial shape, the ROI 
has occurred. The weakness of previously defined quan- 
titative boundaries that separate ROI/non-ROI models 
has lead us to ask if such boundaries exist and, if so, how 
they may be better described. 

We are looking for quantities that distinguish between 
initial conditions that lead to ROI and non-ROI situa- 
tions. It has been suggested that the global quantity A = 
2Tr/rtan = 2 < > / < a?^„ > is a good ROI diagnos- 
tic ( e. g.. iPolvachenko fc ShukhinanI Il981l: iBarnes et al\ 



119861: iMerritt fc Aguilail Il985 t iBellovarv al\ I2008D. 
From analyt ical considerations, iFridman fc Polvachenkd 
(fl98l and IPolvachenko fc Shiikhma^ (IT98ll ) suggest 
that an initial value of Aq > 1.7 will ensure that a sys- 
tem will undergo the ROI. However, the value of A asso- 
ciated with the ROI has bee n foun d to range from » 1.4 
to « 3.0, and lUdrvl (|1993f ) and iHozumi et all (|1996f ) 
discuss situations in which the R OI occurs despite A 
having relatively small values. As iBarnes et al\ (|1986l ) 
speculate, this range is at least partially caused by dif- 
ferent initial distribu t ions o f mass and velocity. Also, 
iPalmer fc Papaloizoul ()1987f ) point out that this quan- 
tity does not have a critical value that separates stable 
and unstable systems, but that the growth rate of the in- 
stability may become small enough as A decreases that 
systems are practically stable. We investigate whether or 
not other global measures can provide a clearer demarca- 
tion between stable and unstable systems. In particular, 
we use the initial virial ratio Qq and the initial fraction 
of mass with radially anisotropic velocities fJ-r,o- As we 
discuss below, A and fir are different ways of measur- 
ing roughly the same characteristic of a system, however 
we focus on fi^ as it has well-defined limits and a simple 
interpretation. 

To clarify the notation that we use in this paper, the 
variables Q, A, fir, fJ-i, and fit refer to the value of those 
quantities at arbitrary times during an evolution. We are 
often interested in the initial values of those quantities, 
which we will denote with a subscript '0'; for example, 
the initial value of the virial ratio is Qq. 

As with previous investigations, we find that two seem- 
ingly different classes of initial conditions lead to the 
ROI; systems with initial velocity anisotropy and dy- 
namically cold, but isotropic systems. The subsequent 
discussions of these two classes focus on their common 
aspects as a way to understand the basic physics of the 



ROI. We will discuss the development of the ROI for the 
two cases in SlU but for the moment we focus on whether 
or not the two cases lead to different global behaviors. 

3.1. Initially Anisotropic Systems 

Models that have enough radial anisotropy undergo the 
ROI; i.e., the shape of the model becomes non-spherical. 
We quantify the presence and strength of the ROI in 
terms of the values of the systemic axis ratios. Our goal is 
to find a link between the amount of radially anisotropic 
mass fir.o initially present in the models and the strength 
of the instability. For each of the initial density profiles 
discussed in § [2l we follow models with Qq = 1.0 (hot), 
0.5 (warm), 0.2 (warm), and 0.1 (cold). Twelve mod- 
els with 0.1 < Ta < 1.2 are created for each Qq, cov- 
ering a spectrum of ^r,o values from nearly isotropic to 
completely radially anisotropic. We have found that the 
radially anisotropic mass fraction has a monotonic, one- 
to-one relationship to the global anisotropy parameter 
A. The exact correspondence depends on the particular 
density profile of the model, but it is not surprising that 
nearly isotropic models with small fir have A ^ 1, while 
models with /i^ ~ 1 have A 3> 1. Thus, the span of /i^.o 
values for our models implies a corresponding range of 
coverage in terms of A. 

We use the innermost 80%-mass axis ratios found in 
each simulation as measures of the strength of the in- 
stability. Specifically, the time at which ^ is minimum 
is found and then the corresponding ^ value is saved. 
While it is common for - and - to have their minimum 
values at the same time, this is not always the case. In 
general, the time of minimum ^ is approximately the 
virialization time (when Q begins to maintain a steady 
value of 0.5). 

Hot and warm models (Qq > 0.2) generally have the 
same behavior, independent of initial density profile (Fig- 
ure [2]). In this figure, the open symbols are derived from 
N ^ 10-* NB0DY2 simulations and the filled symbols 
result from N = 10^ Gadget-2 simulations. There is no 
significant difference between the two. For /ir,o < 0.60, 
the systems remain nearly spherical. Higher amounts 
of initial radially anisotropic orbits lead to non-spherical 
systems. We refer to the firfl values that separate spher- 
ical and non-spherical models as threshold values. For 
models with 0.60 < fir.o < 0.80, the systems take on pro- 
late spheroidal shapes when the instability has maximum 
strength. Systems that are nearly completely dominated 
by radially anisotropic orbits tend to take on fully triax- 
ial shapes at the peak of the instability. The ^ values 
of these heavily radial models tend to decrease as fir,o 
increases, but ^ does not, leading to triaxiality. Over- 
all, there is a slight trend for our models, at all fir.o, to 
become less spherical as they become colder. 

From Figure [H we see that the coldest initially 
anisotropic models with Qq = 0.1 depart from the uni- 
fied behavior of the warm models. As in Figure [2l 
the open and filled symbols are derived from NB0DY2 
and Gadget-2, respectively. The initially constant den- 
sity models barely remain spherical even when mostly 
isotropic, and there is no sudden transformation to sub- 
stantially non-spherical shapes. Overall, ^ values do not 
change much with fir.o, and ^ decreases only slowly as 
fir.o increases. Models with fir.o ^ 0.60 become triaxial 
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when the ROI is at peak strength. The cuspy and Gaus- 
sian models continue to show the previously noted tran- 
sition, but the details are different. First, these models 
are less spherical than their warmer counterparts. Sec- 
ond, the transition to non-spherical systems occurs at a 
smaller value of fir.o ~ 0.40. Third, all systems with 
fir.o ^ 0.40 form nearly prolate systems; there is no 
strong division between prolate spheroidal and triaxial 
systems. Overall, the distinction between models with 
small and large firfl is diminished compared to the warm 
models. This is unsurprising since particles in these mod- 
els have small velocities, and as a result, the difference 
between isotropic and radially anisotropic systems is rel- 
atively minor. The evolution of these cold systems is 
dominated by collapse. 

We have evolved several of these models for much 
longer times to determine the axis ratio behavior on two- 
body relaxation timescales. While some of our mod- 
els reach and maintain their minimum -, - values for 
up to tens of crossing times, all of our models appear 
to become more spherical over a relaxation time. This 
is consistent with the b ehavior of models presented in 
iTheis fc SpurzemI (|1999( l. Those authors speculate that 
the long-term transformation to sphericity stems from 
a violent relaxation process that follows the ROI and a 
subsequent two-body relaxation that occurs between the 
centrally collapsed core (treated as a single very massive 
particle) and particles near the edge of the system. While 
our models do not show a clear distinction between the 
ROI and violent relaxation phases, there is some evidence 
for a weakening of the instability as violent relaxation 
comes to an end. 

From the correspond ence between and ^r,Oi we con- 
firm the speculation in iBarnes et al\ ()1986l ) that some of 
the variation in the value of Aq that separates stable 
and unstable models is due simply to different mass dis- 
tributions. In our models, non-sphericity sets in when: 
^0 ^ 2.1 for initially Gaussian models, Aq > 2.5 for 
initially cuspy models, and Ao > 3.0 for initially con- 
stant density models. We interpret this trend as fol- 
lows. Models with different density profiles have differ- 
ent numbers of particles in their central regions. Any 
bar-like structure that forms near the center of a system 
through the ROI can grow only by ensnaring neighbor- 
ing orbits. In less dense systems, there are fewer parti- 
cles in the central regions, and the bar-like structure will 
thrive only if particles from the outer regions approach 
the center. This occurs more frequently when more of 
the orbits in the outer regions are radially anisotropic, or 
equivalently have smaller /Xr.o(l arger Ap) values. O ver- 
all, our results sup p ort th ose of lBarnes et al\ l|1986( l and 
iMerritt fc Aguilan ()1985f) : larger initial radial velocity 
anisotropy near the center of a system makes it more 
susceptible to the ROI. 

3.2. Cold, Initially Isotropic Systems 

We now turn to the other class of initial conditions 
that lead to the ROI; systems with isotropic velocity dis- 
tributions but extremely small virial ratios, Qo ^ 0.1. In 
particular, we are interested in the behavior of during 
the evolution of systems in the small-Qo cases. Initially, 
Mr,o = 0, but /ir later increases dramatically as any very 
cold collapse leads to 'shell-crossing'. Sets of particles 
with similar initial radial positions fall to the center of 



the potential well and appear to rebound radially out- 
ward where they coexist with other sets of particles that 
are infalling. This overlap of velocity behavior gives rise 
to large radial velocity dispersions which, in turn, in- 
creases f3 in the overlap region . 

Foll owing the work of lUdrvl l|1993f l and lHozumi et all 
()1996f ). we investigate how a) the maximum fi^ achieved 
by the small-Qg cases links to the triaxiality of the re- 
sulting system and b) these systems' axis ratio behavior 
compares to those of initially anisotropic systems. The 
hypothesis is that the maximum /i^ will determine global 
triaxiality and that the /^r-axis ratio relation observed 
in the initially anisotropic systems will continue to hold. 
Five initially isotropic models with Qq = 0.00, 0.02, 0.04, 
0.06; and 0.08 have been evolved for each of the three ini- 
tial density profiles discussed earlier (see § 12. ip . 

3.2.1. pinitial(»") = po 

The initially constant density models undergo such 
strong collapses that a sizeable fraction of particles (~ 
35%) escape the system and leave a mostly spherical 
core for all values of Qq. Except in the coldest systems 
{Qq = 0.00 and 0.02), the axis ratios for these models 
typically remain > 0.85. The coldest models briefly be- 
comes oblate spheroidal with ^ w 0.75. The maximum 
values of fir (/ir,max) for thcsc models appear to be ap- 
proximately 0.5. However, this is a result of coarse time 
resolution of these simultions; fir can reach larger val- 
ues, but only for short intervals < 0.1 tcross(see ii l4.2.1|) . 
Fully isotropic models with slightly "warmer" initial con- 
ditions (Qo = 0.1, 0.2) have similar maximum values of 
fir, mass loss fractions, and axis ratio evolutions. How do 
these models compare to their most similar anisotropic 
counterparts? The Qo = 0.1 models with the weakest 
initial anisotropies {fir.o = 0.30,0.50) have smaller mass 
loss fractions (« 15 — 20%), but still result in a spher- 
ical core. We conclude that the ROI is suppressed in 
these systems simply because sufficient mass occupying 
radially anisotropic orbits is not retained during the col- 
lapse. We point out that not all systems that lose mass 
necessarily remain spherical. If fir,o ^ 0.50, then a sys- 
tem may lose mass and still have enough residual radial 
anisotropy to bring on the ROI. 

3.2.2. pinitiai(r-) oc 

Models with this central density cusp lose approxi- 
mately 15% of their mass after these cold collapses. Un- 
like the constant density cases above, these models all 
reach maximum fir values > 0.85 during their evolutions. 
According to our hypothesis, these systems should all 
undergo the ROI and become triaxial. This is observed, 
but the Qo = 0.08 model remains nearly spherical for 
the first 4 crossing times. After that, the inner 80% of 
the mass becomes prolate while mass exterior remains 
spherical. This change docs not coincide with the mass 
loss which occurs after « 8 crossing times. For models 
with Qo < 0.08, the trends seen in initially anisotropic 
systems, a) weaker anisotropy produces more prolate ob- 
jects and b) stronger anisotropy produces smaller ^ val- 
ues, are also reproduced. The increased central potential 
due to the cusp allows the system to better retain its 
mass through the collapse, which in turn allows larger 
radial anisotropies to develop and lead to the ROI. 
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3.2.3. pinitiai(r-) oc e 

Initially Gaussian density profiles provide a centrally 
concentrated system without a cusp. With our choice of 
Gaussian scalelength (see § I2.ip . the Qo < 0.1 models 
lose < 1% of their mass after their collapses. Like the 
cuspy models above, these systems all reach maximum 
fir values > 0.80 and do not remain spherical. However, 
while the models with Qo < 0.08 follow the previously 
discussed trends (prolate-weaker anisotropy, smaller ^- 
stronger anisotropy), the Qo — 0.08 model remains 
nearly spherical throughout its evolution (^, f > 0.85). 
We have tested the robustness of this result by evolving 5 
different realizations of the Qo = 0.08 model. The max- 
imum variation between axis ratios in the realizations is 
Si 0.1; all realizations stay nearly spherical. 

3.3. Synthesis of Anisotropic and Cold, Isotropic Global 

Results 

We have found that Qo can be a u seful ROI discrim- 
inator for initially isotropic systems. iMerritt fc Aguilarl 
(|l985f ) conclude that Qo < 0.1 should separate ROI sta- 
ble and unstable models that form from isotropic, spher- 
ical collapses. Our boundary value appears to depend on 
the initial density profile, but is « 0.1 for initially cuspy 
(p cx r~^) and Gaussian profiles. For initially constant 
density systems, these cold collapses are apparently too 
violent to allow the ROI to fully develop as large frac- 
tions of mass escape. However, the coldest systems de- 
velop non-spherical shapes before substantial mass loss 
occurs and then return to nearly spherical afterwards. 
We have also evolved systems with larger Qo values. We 
find that during these evolutions, axis ratios typically 
decrease by less than 5% for Qo > 0.1. For Gaussian 
initial densities, the axis ratios decrease by less than 2% 
for Qo > 0.1. Hence, there is a dramatic difference in 
evolution between systems with Qo > 0.1 and Qo < 0.1. 
In this class of models, the global maximum value of 
/ir does not change significantly with Qo, but there is a 
trend for colder collapses to result in larger /Ltr,max- This 
is not surprising and we conclude that the value of /ir, max 
is not a useful diagnostic for the ROI. 

Qo also plays a role in the behavior of the ROI for ini- 
tially anisotropic systems, along with Hr.o- Systems with 
Qo > 0.1 reach similar levels of non-sphericity when they 
have similar amounts of mass on radially anisotropic or- 
bits, independent of initial density profile. In general, 
when firfl ^ 0.60, systems remain spherical. Higher val- 
ues of fj^rfi lead to evolutions that produce prolate and 
triaxial systems. As with the initially isotropic cases, 
Qo ~ 0.1 is again a boundary for ROI behavior. Initially 
anisotropic models with Qo ~ 0.1 have substantially 
lower thresholds to achieve non-sphericity; ^rfi ^ 0.40 
for cuspy and Gaussian models. Constant density mod- 
els remain roughly spherical for /i^.o ^ 0.60 and become 
triaxial with more radial velocity distributions. 

A visual summary of these results is presented in Fig- 
ure m Panels a, b, and c correspond to models with 
constant, cuspy, and Gaussian initial density profiles, re- 
spectively. The shaded areas represent the initial condi- 
tions that give rise to various system shapes at the peak 
of the ROI. For small Qo, every model investigated be- 
came non-spherical. 



4. ONSET AND EARLY DEVELOPMENT OF THE ROI 

We now compare the development of the ROI in 
isotropic and anisotropic systems, describe how the on- 
set of the ROI depends on Qo, and investigate how the 
velocity distributions in these systems change during evo- 
lution. Section [4l1 is intended to complement the discus- 
sion in § 13.11 while Sections 14.21 and [3.21 are counterparts. 

4.1. Initially Anisotropic Systems 

To complement our discussion of the global behavior of 
the ROI for initially anisotropic systems and compare to 
isotropic cases, we now turn to describing the onset of the 
instability in anisotropic models. We have run high time 
resolution versions of a subset of the anisotropic models 
discussed earlier. These evolutions cover approximately 
the first two initial crossing times for each model. We 
again focus on models with Qo-values of 1.0, 0.5, 0.2, 
and 0.1. To cover the same range of initial anisotropy 
as discussed in § 13.11 , models for each Qo were evolved 
for ra = 0.9 (most isotropic), 0.5, and 0.1 (most radially 
anisotropic) . Variations in initial density profiles produce 
qualitatively similar results, which we discuss first. Later 
we compare the behavior of our mode ls to the analytical 
predictions from lHozumi et all ()1996[ ). 

4.1.1. Mass Fraction & Axisratio Behaviors 

Models with = 0.9 have the smallest amounts of 
initial radial anisotropy for each density profile. Initially 
constant density models have fir,o = 0.50, cuspy models 
have fir,o = 0.40, and models with initially Gaussian pro- 
files have fj,r,o = 0.10. For all values of Qo, the value of 
fir rapidly decreases up to the time of maximum collapse 
(when the system is in its most compact configuration), 
with fit increasing to a maximum during the same pe- 
riod. After maximum collapse, fir values rise again with 
colder models showing stronger and more rapid increases. 
In the coldest systems (Qo < 0.2), fir can reach and sur- 
pass the threshold values discussed in § 13.11 However, 
the value of the centrally-concentrated, isotropic mass 
fraction fii^c remains greater than about 0.20. During all 
of these evolutions, the axis ratios maintain quite spher- 
ical values (| > 0.90 and ^ > 0.85). Deviations from 
sphericity correspond to times during which fir is near 
or above its threshold value and fii c is near its minimum 
value. 

More initial radial anisotropy occurs in systems with 
ra = 0.5; fir,o=0-90, 0.80, and 0.60 for models with ini- 
tially constant, cuspy, and Gaussian density profiles, re- 
spectively. As with the more isotropic models, the fir 
values drop to their minimum levels just before maxi- 
mum collapse in these evolutions. The fii and fit val- 
ues both increase during this time, with warmer mod- 
els showing larger fit levels and colder models having 
larger fii values. After maximum collapse, fir values in- 
crease to at least threshold levels for all Qo values. In 
the warmer models, fii^c values remain larger than 0.20 
during the interval when fir is above the threshold. Axis- 
ratios for these warm models do not change significantly 
from spherical values (« 0.90) during this part of the evo- 
lution. Colder models have fii c ^ 0.20 while fir is above 
the threshold, and their axis ratios decrease while both 
of these conditions hold. Overall, the amount of axis ra- 
tio change increases with decreasing Qo even as the time 
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at which the axis ratios reach their minimum values de- 
creases. Minimum axis ratios are typicaUy reached after 
several maximum collapse times have passed. Constant 
density models find their axis ratio minima earlier than 
cuspy models. Axisratios in Gaussian density models 
take the longest amount of time to reach their minimum 
values, around 20 maximum collapse times. 

When Ta = 0.1, the models are initially completely 
radially anisotropic for all density profiles. Since no 
mass is isotropic, our previous findings would suggest 
that these systems should immediately begin to change 
shape. However, no models show substantial changes in 
axis ratios until after their times of maximum collapse. 
Minimum axis ratio values are typically reached between 
4 and 8 maximum collapse times, independent of ini- 
tial density profile. In these extreme models, the time 
evolution of for initially constant and cuspy systems 
behaves somewhat differently than for Gaussian systems. 
In constant and cuspy models, decreases to at least 
threshold values by the time of maximum collapse and 
then increases again. The /i^ values in the Gaussian mod- 
els remain above threshold values for at least the first two 
crossing times. The outer edges of constant and cuspy 
models develop significant levels of isotropy (/i^ > 0.40), 
while Gaussian models do not show this behavior. De- 
spite these differences, all models share the following be- 
havior; the amount of centrally-concentrated, isotropic 
orbits increases during an evolution, and the increase 
is greater for colder initial conditions. In general, ^i^c 
reaches a value of 0.20 between 1.5 and 2.0 tcross- Again, 
as discussed previously, the strength of the axis ratio 
change in a given system increases with decreasing Qq. 

4.1.2. Velocity Dispersion Evolutions 

iHozurni et al\ (119 961 present predictions [based on 
work in lKan-va eF* i l. (1996.)] for the evolution of veloc- 
ity dispersions in systems with initially i) isotropic ve- 
locity distributions and ii) constant and cuspy (p oc r~^) 
density profiles as well as results of numerical models. 
We directly compare our results to theirs for analogous 
models and note resemblances for our Gaussian models, 
which do not have analytical counterparts in their work. 

Tangential dispersions (atan) in our initially constant 
density and cuspy = 0.9 models (which are most sim- 
ilar t o the completely isotropic models in iHozumi et al\ 
|1996[ ) evolve at least as fast as radial dispersions (ur), 
in qualitative agreement with analytical models. In par- 
ticular, colder (Qo = 0.1 and 0.2) models have disper- 
si on evolutions similar to their numerical counterparts 
in lHozumi et al\ ()1996f ). Our Gaussian models have dis- 
persion evolutions that are qualitatively similar to those 
for the cuspy models. For example, Figure [5] shows how 
the dispersions in a cold (Qo = 0.1), cuspy model evolve 
up to the time of maximum collapse and is quite similar 
to Figure 5b in lHozumi et al\ (|1996( ). We note that our 
models disagree with the analytical predi ctions in much 
the sa me way that the numerical models in lHozumi et al\ 
(|1996[ ) do. The lack of inifinite contraction in numeri- 
cal models forces them to behave as though they have 
roughly constant density cores, which analytical models 
predict should evolve isothermally. The analytical mod- 
els also neglect system edges, which leads to numerical 
models with smaller than predicted dispersions in their 
outer regions. While we find good agreement with the 



previous numerical models, there are two aspects of our 
evolutions that are different. 

First, our warm models (Qo = 1-0 and 0.5) do not 
maintain small isotropic cores as maximum collapse is 
achieved. Instead, these systems become almost com- 
pletely ta ngentially anisotropic by the time of maximum 
collapse. iHozumi et al\ (|1996f ) only considered models 
with Qo ^ 0.3, so this behavior has not been described 
before. Second, for all our models, the decrease in radial 
dispersions in the outer regions is more severe, in magni- 
tude and abruptness, than the previous models. Rather 
than a smooth transition between inner regions where Ur 
increases to outer regions where (t^ decreases, we see a 
"break" radius that divides the inner and outer regions. 
This break radius moves inwards as the time of maximum 
collapse is approached. 

The break radius arises from the initial conditions we 
have set in our models. Radial dispersions are created by 
giving particles both positive and negative radial veloci- 
ties. Particles near the edges of systems have the largest 
anisotropics, and so, radial velocities as well. Particles 
that have positive radial velocities quickly expand the 
system. When mass shells are then created, the outer- 
most shells arc comprised entirely of expanding parti- 
cles; the average velocity is large and positive, but the 
dispersion is nearly zero. As time passes, particles with 
positive radial velocities from inner regions expand to a 
point where they are surrounded only by particles with 
similar motions, decreasing the dispersion. In this way, 
the outer low-radial dispersion region grows inward with 
time. By the maximum collapse time, the break radius 
is replaced by a smooth transition between the inner, Ur- 
increasing, regions and outer, dr-decreasing, regions. 

Models with more radially anisotropic velocity distri- 
butions, i.e., those with — 0.5 and = 0.1, also 
behave qualitatively the same as their more isotropic 
cousins. Tangential dispersions evolve very similarly 
to the predictions of analytical models; they increase 
throughout the system, but the largest increase occurs 
at the center. Radial dispersions increase at the centers 
of systems to maintain near-isotropy and decrease near 
the edges. The biggest difference between these models 
and the = 0.9 versions is the magnitude of the dis- 
persions. For the cuspy model with Ta = 0.1, the initial 
radial and tangential dispersion values are roughly con- 
stant in radius, but cr^ ~ SOdtan- So. even as crtan grows, 
the systems can stay quite radially anisotropic. 

4.1.3. Unified Picture of ROI Onset in Initially Ansotropic 

Systems 

A common story emerges from comparing models with 
different initial density profiles and initially anisotropic 
velocity distributions. The amount of centrally- 
concentrated mass on isotropic orbits plays an important 
role, as previously hypothesized in iMerritt fc Aguilarl 
(|198S ). However, the /i^.c value, by itself, does not corre- 
late with a system undergoing the ROI. We have found 
that a combination of and fii^c can explain our simu- 
lations in the following way. If iii^c ^ 0.20 and reaches 
levels larger than the appropriate threshold value, then 
the axis ratios decrease rapidly. We stress that, individ- 
ually, the two conditions do not describe the onset of the 
ROI. Larger values of give rise to larger changes in 
axis ratios, as long as ^i^c is small enough. The = 0.1 
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models give an interesting caveat to our picture. These 
models begin totally radially anisotropic, but the sub- 
stantial changes in axis ratios do not begin until around 
the time of maximum collapse. 

Keeping with the picture of the rad i al orb it instabil- 
ity suggested by iPalmer fc Papaloizoul ()1987t ) , any non- 
spherical distribution of mass creates torques on orbits 
not aligned with the major axis of t he system. We note 
that this idea also has precedents inlLvnden-Belll (119791 ) 
and iFridman fc Polvachenkcl ()1984f K These torques 
would be largest when the system is most compact. In 
Appendix [X\ we discuss a model of this situation; a test 
mass and two off-center masses. This toy model results 
in a torque on the test mass that increases as the test 
mass radius decreases, which we contend will hold for 
more realistic bars as well. As a result of torque increas- 
ing during the collapse, more orbits will be trapped by 
the bar and further strengthen the bar. 

This idea helps explain why the warm = 0.1 models 
take longer to form weaker bars than the cold models 
with the same radial anisotropy. In the warm models, 
only a fraction of the mass will converge to the center, so 
the torques that result are smaller and require a longer 
time to have an effect. Colder models allow more mass to 
collapse, providing a stronger bar and torques. We note 
that this picture also qualitatively explains why more 
isotropic models form weaker bars. When the orbits of 
particles are not nearly radial, and have appreciable an- 
gular momentum, the torques exerted cannot align those 
orbits to the bar as quickly. 

The question becomes, why don't these systems be- 
come extreme bar s? Our pictu r e sup ports the hypoth- 
esis forwarded bv lAdams et al\ (|2007| ). namely that an 
instability of individual orbits along principal axes in 
triaxial systems isotropizcs the orbits and halts bar 
growth. In looking at the anisotropic mass fractions, 
we have pointed out the important role of the centrally- 
concentrated isotropic mass, iii^c- The value of this quan- 
tity is closely linked to the growth of a bar. In cases 
where ^ is sufficiently low and is high, bar growth 
occurs, followed by an increase in c- We equate this 
growth in fXi^c with the outcome of the orbit instabilities 
present in the triaxial system. 

4.2. Cold, Initially Isotropic Systems 

The hypothesis [also discussed by iHozumi et al\ 
(|1996f )] is that fir for an initially isotropic system reaches 
the level indicated by the initial anisotropic models at 
the same time the system begins to deviate from spher- 
ical symmetry. In effect, the isotropic system develops 
radial anisotropy through collapse beyond which point it 
behaves like an initially anisotropic system. To test this 
idea, we look at high time resolution versions of the cold, 
isotropi c models discu ss ed in § 13.21 

Like iHozumi et ai\ (|1996l ). we find that initially 
isotropic systems with Qq < 0.1 remain nearly, but not 
perfectly, spherical during their collapse phases. Specif- 
ically, the axis ratios slowly decline from their initial 
values until just before the point of maximum collapse, 
when the half-mass radius is minimum. We attribute 
this slow initial change in axis ratio values t o the type of 
evolution discussed in iLin. Mestel. fc Shul ()1965f ). Any 
slight departure from sphericity becomes magnified dur- 
ing gravitational collapse. The axis ratio change is ex- 



pected to increase as the rate of collapse increases, or 
alternatively, as the system gets colder. This is indeed 
what is seen for all of our models, independent of density 
profile, but the effect is more prominent in the initially 
constant density models. As an example of this trend 
from the constant dentisy models, the Qq = 0.08 model 
has A- = —0.04 during the time until maximum col- 
lapse, while the Qq = 0.00 model has = —0.18 over 
the same amount of time. The ^ values given here refer 
to the axis ratios of the innermost 95% of the mass. The 
80%-mass axis ratios show similar, but more exagger- 
ated, behavior since their initial axis ratio values tend to 
be initially farther from 1.0, compared to their 95%-mass 
counterparts. We now continue with more detailed de- 
scriptions of the early evolutions of these systems (from 

i = to t = 2 icross)- 

4.2.1. Pinitial(r) = po 

In general, the constant density cases maintain a large 
degree of isotropy throughout their infall periods. We 
consider the infall period to be the time from the be- 
ginning of the simulation to the maximum collapse time, 
tnic ~ 1-1 Across for thesc models. The Qq = 0.00 and 
Qo = 0.02 models maintain isotropy in at least the in- 
ner 50% of the mass in the system continuously until the 
maximum collapse time. For models with 0.04 < Qq < 
0.08; fii > 0.50 until « O.ltcross before the maximum col- 
lapse time. Between this time and the maximum collapse 
time, the central isotropic mass fraction c decreases as 
more and more orbits become tangentially anisotropic. 

At the maximum collapse time, all systems rapidly de- 
velop large amounts of radial anisotropy, with a corre- 
sponding decrease in central isotropy. This is due to 
particles rebounding outward from their collapse occu- 
pying the same volume as other particles which are still 
infalling. This period of radially anisotropic domination 
is very short-lived. The duration is about 0.1 icross for 
the Qo = 0.08 model, and decreases with Qo- After the 
brief period of radial anisotropy, /ii_c increases to its pre- 
maximum collapse value. 

When /ii.c briefly decreases, the axis ratios change 
abruptly (for both the 80%-mass and 95%-mass ver- 
sions). These changes are not monotonic; the axis ratios 
can increase and decrease by sizeable amounts, but by 
the time fir reaches its maximum, the axis ratios have de- 
creased to values lower than they were at maximum col- 
lapse time. This overall decrease scales inversely with the 
Qo-value. For the coldest models, axis ratios reach their 
minimum values soon after the maximum collapse time. 
Once iii c returns to its pre-maximum collapse value, the 
axis ratios do not change appreciably until mass begins 
to escape from the system. This occurs fairly rapidly for 
the Qo = 0.00 model (mass loss begins at w 1.6tcross), 
but takes longer to develop (after 2tcross) for models with 
larger Qq. 

We find that the velocity dispersions in the models 
with Qo > 0.02 evol ve similarly to the c onstant den- 
sity model discussed in lHozumi et all (jl996l . compare our 
Figure [6] with their Figure 3a) . Initially, both the ra- 
dial and tangential dispersions are equal and constant 
throughout the system. As the system evolves to about 
2/3 of the maximum collapse time, the dispersions in the 
central regions grow contemporaneously, while the radial 
dispersion grows much more slowly than the tangential 
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dispersions beyond about the half-mass radius. During 
the final third of the evolution to maximum collapse, 
the innermost dispersions grow much more rapidly than 
those farther out in the system. The central half of the 
system is isotropic while the outer half is very tangen- 
tially anisotropic. The Qo = 0.00 model on the other 
hand shows a diffe rent behavior, i nuch more like the an- 
al ytical results in iHozumi et al\ (|1996D [based on work 
in iKan-va et al\ ()1996[ )]. The radial and tangential dis- 
persions evolve together for approximately the first 3/4 
of the maximum collapse time, maintaining system-wide 
isotropy. As the system continues towards maximum col- 
lapse, the radial dispersion begins to lag behind the tan- 
gential dispersions; approximately 30% of the outermost 
mass is on tangcntially anisotropic orbits. 

4.2.2. pinitiai(r-) oc r"^ 

These models maintain central isotropic regions for ap- 
proximately 0.5 tcross- The maximum collapse times for 
these models are imc ~ 1-1 Across- For the Qo = 0.00 case, 
the decrease in central isotropy after 0.5 tcross is accom- 
panied by a growth in radial anisotropy with very little 
tangential anisotropy present. Models with Qo > 0.02 
also lose their central, isotropic regions, but they develop 
a substantial amount of mass on tangcntially anistropic 
orbits (/it = 0.8). The loss of central isotropy persists 
for 0.5 tcross ^ t ^ 1-5 tcross for the Qo = 0.00 and 
Qo — 0.02 models. As Qo increases, it takes longer for 
the central isotropy to return to pre-maximum collapse 
levels. In the Qo ~ 0.08 model, ^i^c drops to and re- 
mains around 0.15 until about 4 tcross- For those models 
that develop large amounts of tangential anisotropy, it 
disappears prior to maximum collapse. As the tangcn- 
tially anisotropic mass fraction decreases, /i^increases, 
typically reaching its maximum value shortly after max- 
imum collapse. 

The axis ratios (both 80%-mass and 95%-mass) for 
these systems show slow decreases (like those expected 
simply from contraction) until the /i^. values reach the 
threshold level discussed in § 13.11 At this point, fii^c is 
at its minimum and the axis ratios decline more sharply. 
This decline continues until fii^c ^ 0.25. As discussed 
earlier, this minimum in the central isotropic mass frac- 
tion lasts for shorter periods of time for models with 
smaller Qo. At the same time, the total axis ratio change 
is inversely proportional to Qo, so the axis ratio decrease 
is most rapid for lower-Qo systems. Overall, the axis ra- 
tios in the coldest systems reach their minimum values 
earlier (« 2 tcross) than those in warmer systems (after 

4 tcross)- 

A longer simulation of the Qo = 0.08 model with the 
same time resolution gives some further insight. The 
95%-mass axis ratios remain at nearly spherical values 
for the first 8 crossing times (until mass loss occurs). 
The 80%-mass axis ratios remain nearly spherical until 
~ 4 tcross, after which they decline together indicating 
a moderately strong bar has formed. The axis ratios 
(both 95%-mass and 80%-mass) for the initially isotropic 
Qo = 0.1 model remain above 0.9 until a mass loss event 
occurs. We have evolved an additional cuspy model with 
Qo = 0.09 to better pin down the boundary of systems 
that are susceptible to the ROI. The axis ratios for the 
Qo = 0.09 model behave similarly to the Qo = 0.08 case. 
We interpret this behavior as evidence that Qo = 0.1 



is very nearly the value that separates ROI stable and 
unstable systems for this density profile. When Qo is 
near 0.00, the instability occurs simultaneously through- 
out a system. Systems with borderline-unstable Qo val- 
ues show delayed onset of the instability and once the 
instability begins, only the inner part of the system ap- 
pears to take part. 

As with the constant density models, the evolution 
of the velocity dispersions in the cuspy systems with 
Qo ^ 0.02 a lso agree with the beh avior of similar systems 
discussed in lHozumi et al\ ()1996[ ). Until about half of the 
maximum collapse time, the systems roughly follow this 
behavior: dispersions in the central regions grow rapidly 
but isotropically, the tangential dispersions increase more 
slowly near the edges, and the radial dispersions remain 
nearly constant outside the central « 10% of the mass. 
The Qo = 0.00 model differs during these early stages 
of the collapse in that the dispersions evolve isotropi- 
cally throughout the system. As the systems continue 
to maximum collapse, radial anisotropy begins to dom- 
inate from the center outward (independent of Qo)- By 
three-quarters of the maximum collapse time, the cen- 
tral 30-50% of the mass is on radially anisotropic or- 
bits. At the maximum collapse time, the central-most 
5-10% of the mass is i sotropic, with t he rem ainder being 
radially anisotropic. iHozumi et al\ ()1996f ) explain the 
isotropic evolution of the central dispersions by stating 
that the density remains nearly constant (spatially) in 
these regions; isotropy is then the expected behavior (see 
§ 14.2. ip . We find that the density profile slowly changes 
from the initial r^^ shape to a near-de Vaucouleurs shape 
at the maximum collapse time, maintaining a cuspy cen- 
ter (even if it is not exactly r~^). 

4.2.3. Pinitial('-) OC e"''^ 

Compared to the previous systems, these models have 
the shortest maximum collapse times, t^c ~ 0.5 tcross - 
The central isotropic mass fraction begins to decrease 
about 0.1 tcross prior to the maximum collapse. Like 
the situation for cuspy profiles, the Qo — 0.00 model 
mainly develops strong radially anisotropic motions as 
the isotropy is lost; only a small fraction of mass (« 10%) 
is tangcntially anisotropic. Models with Qo > 0.02 lose 
central isotropy as well, but they rapidly become near- 
completely tangcntially anisotropic before maximum col- 
lapse. This tangential anisotropy is quickly lost as the 
system develops large radially anisotropic mass fractions 
(« 0.50) by the maximum collapse time. The fir val- 
ues reach maximum values of around 0.80, independent 
of Qo, approximately 0.1 tcross after maximum collapse. 
The central isotropic mass fraction dips to 0.05 before 
maximum collapse, typically when most of the mass is 
on tangcntially anisotropic orbits. For the Qo = 0.00 
case, the value remains at this minimum value during 
the entire rise of fir- In models with Qo > 0.02, fii^c in- 
creases slightly by the time of maximum collapse. The 
amount of increase is proportional to Qo; for example, 
fii c increases to 0.10 for Qo = 0.02 while fii c increases 
to 0.20 for Qo = 0.08. 

The behavior of axis ratios in these models is similar 
to those for the cuspy models. The models present slow 
initial declines in axis ratios until the radially anisotropic 
mass fraction increases above the threshold value for this 
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density profile 0.5), occurring at very nearly the time 
of maximum collapse. At this point in time, the axis ra- 
tios start to decline more rapidly. This rapid decline con- 
tinues until reaches its maximum and while fii^c is less 
than « 0.20. As before, the total decrease in axis ratio 
values is inversely proportional to Qq. Since all models 
reach /Zr,max at approximately the same time, colder ini- 
tial conditions produce steeper axis ratio changes. Like 
the cuspy models, we again see a trend for the axis ra- 
tios in the coldest models to reach minimum values ear- 
lier (between 1 and 2 tcross) than warm models (after 3 

Across) ■ 

Velocity dispersions for these models evolve differently 
than in either the constant or cuspy density profile cases. 
The Qo = 0.00 model dispersions stay isotropic for 
roughly half the time to maximum collapse. By 75% 
of the maximum collapse time, the radial dispersion has 
a) evolved with the tangential dispersions at the very 
center, b) lagged behind tangential values out to the 
half-mass radius, and c) reached values larger than tan- 
gential beyond the half-mass radius. At the maximum 
collapse time, the system is isotropic at the very center 
and around a radius about twice the half-mass radius, 
but is radially anisotropic everywhere else. Models with 
Qo > 0.02 have different dispersion behaviors. The ra- 
dial dispersions actually decrease in time over most of the 
radial extent of the systems. Since the tangential disper- 
sions increase during the same time, these systems de- 
velop large amounts tangential anisotropy. The amount 
of radial dispersion decrease is linked to Qo; lower Qo 
leads to smaller decreases. The very centers of all of 
these models remain isotropic throughout the maximum 
collapse time. By the maximum collapse time, the radial 
dispersions have grown faster than the tangential disper- 
sions inside the half-mass radius but continue to lag them 
outside this radius. 

4.2.4. Unified Picture of ROI Onset in Cold, Initially 
Isotropic Systems 

As previously discussed for initially anisotropic sys- 
tems (S I4.1.3p . the roles of fir and fii_c arc key to under- 
standing if and how the radial orbit instability will occur. 
If fir rises above the appropriate threshold value and fii^c 
is small (< 0.20), then a system will change its shape 
dramatically. The exact amount of axis ratio decrease 
depends on the difference fXr — fJ-i,c during this interval. 
Colder initial conditions lead to smaller minimum fii^c 
values, and larger axis ratio changes. Then, independent 
of the fir value, when fii^c increases above « 0.20, the 
axis ratios do not substantially change further. 

As discussed in § 13.11 the threshold values of jU^ for 
systems with 0.2 < Qo < 1.0 are: « 0.65 for constant 
density, w 0.60 for p cx r~^, and « 0.55 for Gaussian 
initial density. The initially isotropic systems develop 
these fir values at approximately the same times that 
the axis ratios begin to decrease in every case. For ex- 
ample, the Gaussian density system with Qo = 0.00 
has a well-defined axis ratio drop that starts just be- 
fore t = 0.5 tcross- In this model, /i,. = 0.55 when 
t = 0.48 tcross- The axis ratio and anisotropic mass frac- 
tion behaviors early in the evolution (t < 2 tcross) are 
shown in Figure [71 The vertical dashed lines indicate 
the times at which /i,. = 0.55 (left) and fii^c ~ 0.20 
(right). Note the changes in the axis ratio behaviors 



at these times. In general, the axis ratios for cold, ini- 
tially isotropic models reach their minimum values ear- 
lier (after approximately 2 or 3 maximum collapse times) 
than those for initially anisotropic models (after approx- 
imately 5-10 maximum collapse times). 

5. SUMMARY & CONCLUSIONS 

We have created a variety of self-gravitating, iV-body 
systems to investigate the behavior of the radial orbit 
instability (ROI). Our models use modest particle num- 
bers (N — 10'*, 10^), but have been shown to simulate 
collisionless behavior for at least one two-body relaxation 
time. Tests with various softening length/particle num- 
ber combinations lead us to the conclusion that minor 
variations in these quantities do not play significant roles 
in determining the evolution of our models. With fixed 
N and softening length, we have created a suite of models 
with various a) initial virial ratios Qo (relating the total 
kinetic and potential energies), b) levels of initial veloc- 
ity anisotropy, and c) initial density profiles (constant, 
cuspy, and Gaussian). 

In a global sense, the behavior of the ROI, in terms of 
its strength and overall effect on system shape, is closely 
linked to the fraction of system mass that follows radially 
anisotropic orbits, fir- We find very clear threshold val- 
ues of firfi in our simulations of initially anisotropic sys- 
tems that divides models into spherical and non-spherical 
cases. The exact threshold value of firfi depends on the 
initial density profile of the model, but when ^^,0 ^ 0.50, 
the warm and hot systems investigated remain nearly 
spherical throughout their evolutions (at least 20 cross- 
ing times). With a modest amount of additional radial 
anisotropy, our models transform from spherical shapes 
to very prolate, bar-like shapes. Higher radial anisotropy 
leads to systems that take on triaxial, bar-like forms. We 
find that the /i^.o quantity is a reliable predictor of the 
final shape of the system, as long as the system is suffi- 
ciently warm (Qo > 0.2). 

In terms of initial anisotropy, dynamically cold systems 
behave quite similarly to one another. This is unsurpris- 
ing as the systems' lack of substantial kinetic energy nec- 
essarily relegates any differences in velocity distributions 
to play a minor role. In other words, when particles can 
only have slightly non-zero speeds, it doesn't matter how 
those velocities are oriented. The early evolution of such 
systems will be dominated by gravity, and the dominant 
motion of all particles will be infall. Figure |4] represents 
a summary of our findings regarding the global behavior 
of the ROI, for both initially isotropic and anisotropic 
models. 

Looking at cold, initially isotropic systems, we have 
found that the ROI sets in when a) the amount of mass 
on radially anisotropic orbits is larger than the threshold 
value and b) there is little mass on isotropic orbits at 
the centers of the systems. In these models, any radial 
anisotropy present must develop from free-fall motions. 
During initial collapse the mass acquires radially inward 
velocities. However, since the particles fall together, no 
large dispersions develop. Only after some of the mass 
has "rebounded" away from the center do systems de- 
velop large amounts of radial anisotropy. Accordingly, 
the point of maximum radial anisotropy occurs after the 
time of maximum collapse of the system. 

The onset of the ROI in systems with initially radial 
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anisotropic distributions is also tied to the amount of 
radially anisotropic and centrally-concentrated isotropic 
mass. We again see that the ROI begins when Hr is larger 
than the threshold value and fii^c is less than about 0.20. 
However, even in systems that arc completely radially 
anisotropic to begin with, the ROI will not begin until 
after the system has maximally "collapsed" . This is not a 
true collapse, like in the cold isotropic systems, but refers 
instead to a measure of the average radius of the system 
reaching its minimum value. Typically, some fraction of 
the mass in the system is actually expanding at this time. 

This behavior leads us to suggest that collapses, the 
shrinking of at least some part of a system, play a key 
role in the development of the ROI. When the systems 
are most compact, torques caused by bar-like mass distri- 
bu tions are strongest. Keeping with the idea put forward 
in iPalmer fc Papaloizoul (|l987t ) , these stronger torques 
can align more orbits with the bar more quickly than if 
the system were more extended. The collapses allow the 
torquing mechanism of the ROI to act more efficiently. 

Our models also serve as evidence that the isotropy- 
creating orbit instability described by lAdams et al\ 
()2007t l is involved as the ROI is halted. Note that this 
isotropy-creating instability is completely different from 
the ROI. We do not find extreme bars among our models; 
the smallest ^ value is w 0.4. So, something is halting the 
torquing mechanism discussed earlier. Certainly, a sys- 
tem expansion would weaken the torques produced, just 
as contraction strengthened them. The bar would stop 
growing simply because it could not continue to add more 
orbits to its distribution. However, if this scenario alone 
describes how the ROI stops, there would be no change 
to the velocity properties of particles in the bar after it 
had grown; particles that rolled along the bar previously 
would continue their motions. However, we see that the 
centrally-concentrated isotropic mass fraction increases 
after the bar forms. This behavior is nicely explained by 
the Adams et al. instability. A strong triaxial (or at least 
spheroidal) mass distribution is formed. Any orbit mov- 
ing along a principal axis is unstable to being deflected 
away from its current trajectory, turning radial orbits 
into isotropic ones. These isotropic orbits will make the 



central regions more spherical, further weakening the bar 
and the torques it can cause. 

Figure [S] illustrates one particle's motion along the in- 
termediate axis of the bar (formed by the ROI) as it 
changes after the time of maximum collapse (t sa 1 Across)- 
This particle was chosen to lie close to a plane that runs 
through the center of the system perpendicular to the in- 
termediate axis of the final bar. Additionally, this parti- 
cle has a small velocity in the intermediate axis direction, 
repres enting our closest ap proximation to the orbits stud- 
ied in lAdams et al\ (|2007[ ) . Note that the position and 
velocity were determined after the maximum collapse 
time. The curve shown represe nts the growth factor y de- 
scribed in lAdams et al\ ()2007f ): the particle's intermedi- 
ate axis position throughout the simulation divided by its 
initial position. Exponential growth results in linear seg- 
ments. We see that there is a chaotic, but roughly linear, 
rise starting att = 1 icross followed by a slowly oscillating 
but non-increasing segmen t. This behav i or is s imilar to 
that shown in Figure 4 of lAdams etHU (l2007h . While 
the likeness is rough at best, we remind the reader that 
our system is still undergoing axis ratio cha nges through- 
out th e system during this interval unlike in lAdams et al\ 
(|2007f) where the axis ratios are spatially and temporally 
fixed. This type of behavior occurs in our simulations 
for orbits with a wide range of sizes and orbit timescales. 
Together with the increasing central isotropy discussed 
earlier, the presence of the Adams et al. instability sig- 
nature in our simulations suggests that the ROI can be 
halted by such an instability. However, the number of 
orbits that are vulnerable to the Adams et al. instability 
is relatively small as particles must move nearly along 
principal axes. Whether or not this population of orbits 
is suHicient to actually halt the ROI must be a topic of 
further investigation. 
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APPENDIX 
TOY MODEL OF TORQUES 

We create a simple situation to illustrate our assertion that the torque experienced by a particle with mass much 
smaller than the mass in a bar-like distribution will increase as the overall system shrinks. Our model consists of a 
test mass and two larger masses, placed symmetrically along the x-axis. Figure |9| shows the geometry of the situation. 
For simplicity, we consider a planar orbit only. 

The net force exerted on the test mass rrit by masses mi and m2 is, 

'mirit , m2r2t 
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'2t 



Masses mi and m2 are each located a distance £ from the origin (ri = 7'2 = ^)- so 
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Using the fact that 6i 



92, we can re- write these distances as, 
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Terms are grouped this way because we make the argument that as the system is contraeting (or expanding), the bar 
length should change roughly in the same way as rt- In what follows, we assume that the quantity (./rt will remain 
constant during collapse. 

The torque about an axis through the origin (perpendicular to the page) is just r = r t x i^net- With the net force 
from Equation I AH we have 

■"mirt X Tit , mart x r2t 



f = —Grrit 
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We can simplify the cross-products in this expression by noting that 
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so Equation I A4I simplifies to 



rit = rt - ri 



r2t = n - r2. 
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These cross-product terms point in opposite directions; in this setup, ri x rt is in the —z direction and r2 x rt is in 
the +z direction. We use the relationship between Oi and 62 to write the magnitude of the cross-products as 



n X rt\ 



rirt sin 62 = £rt sin 62 



and 



Taking mi = 7712 = to. Equation 



\r2 X rt \ = r2rt sm6'2 
5] can now be written as 

T = —Grritmirt sin 62 
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Using the expressions in Equation IA3[ we reach our final expression for the torque on the test mass in this geometry. 
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The terms with i/rt or 62 just 



describe the geometry of the situation, but the lone 1/rt determines how the torque on the test mass changes with 
radial position. If the test mass is near the j/-axis {62 ~ '"■/S), then the torque exerted will be small independent of rt 
since /i and /2 will be approximately equal. Test masses near the a;-axis are already in the vicinity of the bar {62 ~ 0), 
and also have very little torque exerted on them. Masses at intermediate angles will be most affected by the torque. 
More realistic models could be constructed by adding more pairs of masses to the system. Superimposing pairs would 
add to the complexity of the geometry term, but each pair would still have a 1/rt dependence. 
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Fig. 1. — We illustrate important initial condition information for an anisotropic system with a cuspy density profile (p oc r~^) and 
Qo = 0.5. The log-log density profile of the system is shown in panel a. Panel b shows the anisotropy profile. The average spherical 
component velocity profiles and velocity dispersion profiles are illustrated in panels c and d, respectively. 
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Fig. 2. — Minimum axis ratios for the innermost 80% of mass in models that have initial radial anisotropy and virial ratios Qo > 0.2. In 
each panel, models with Qo = 1-0 are denoted by square symbols, Qo = 0.5 models are marked with triangles, and systems with Qo = 0.2 
are represented by diamonds. The empty symbols are derived from the NBODY2 simulations with A*^ = 10* particles, while the filled 
symbols represent Gadget-2 simulations with N = 10^ particles. The similar patterns followed by initially constant density models are 
shown in panels a) and b). Panels c) and d) show the similarities between axis ratio behavior in initially cuspy models. The same basic 
patterns also appear for initially Gaussian density models, panels e) and f). As the initial amount of radial velocity anisotropy increases, 
models become non-spherical. 
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Fig. 3. — Minimum 80%-mass axis ratios for cold Qo = 0.1, initially anisotropic systems. As in Figure [J] the open and closed symbols 
refer to evolved with NBODY2 and Gadget-2, respectively. The cuspy and Gaussian models (denoted by triangles and squares, respectively) 
follow the same basic trend for both ^ and ^ . Models with initially constant density (represented by diamonds) do not show much variation 
with /ir, but these systems represent the remnants of models with significant mass loss. 
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Fig. 4. — A visual summary of the global behavior of the ROI for initially isotropic and anisotropic models. Panels a, b, and c correspond 
to models with constant, cuspy, and Gaussian initial density profiles, respectively. Models with initial conditions in the dark gray areas 
will evolve to triaxial shapes at the peak of the ROI. Initial conditions in the light gray region result in prolate systems. Systems that start 
from the unshaded region will remain more-or-less spherical throughout their evolutions. For the constant models in panel a, note that 
the models with Qo ^ 0.1 had significant amounts of mass loss. The shape of the final virialized system refers only to the remnant, and 
drawing strong conclusions about the ROI from these cold models is problematic. 
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Fig. 5. — Comparisons between spherical velocity component dis persions in a cuspy, in itially radially anisotropic model with = 0.9 and 
Qo = 0.1. The overall behavior of the curves matches Figure 5b in lHozumi et ai\ II1996I ). The curves represent dispersions at four different 
times during the evolution up to the time of maximum collapse fmc ■ The earliest set of curves is essentially the beginning of the simulation, 
one-hundredth of the time to maximum collapse. The next set is halfway to collapse i = 0.5 tmc, the third set is for t = 0.75 tmc, and the 
upper-most curves show the situation at the time of maximum collapse. The vertical lines show the location of the half-mass radius at each 
interval. Note that the half-mass radius decreases as the time of maximum collapse is approached. With this level of initial anisotropy, the 
axis ratios do not reach their minimum values until many maximum collapse times have passed. 
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Fig. 6. — Comparisons between spherical velocity component dispersions in the initially isotropic model with Qo = 0.04. The overall 
behavior of the curves matches Figure 3b in lHozumi et al\ iwM) . The curves shown have the same meaning as those in Figure [5] Again, 
take note that the half-mass radius marks move to the left in the figure as time increases. Axisratios typically reach their minimum values 
after roughly two — three maximum collapse time for these models. 
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Fig. 7. — The top panel shows the axis ratios as a function of time in the initially isotropic, Gaussian density profile system with Q = 0.00. 
The behavior of the velocity mass fractions is plotted in the bottom panel. The vertical dashed lines mark the times at which a) left line, 
55% of the system mass is on radially anisotropic orbits and b) right line, 20% of the central mass follows isotropic orbits. 
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Fig. 8. — The logarithm of the growth factor y for an orbit in a high-time resolution simulation with Qo = 0.00 and an initially cuspy 
density profile. This particle is moving in a plane perpendicular to the intermediate axis of the bar at t = 1 tcross (when the growth factor is 
1). The orbit then experiences a roughly exponential growth perpendicular to the bar after which it settles into a quasi-periodic oscillatory 
phase. This behavior is reminiscent of orbits that undergo the instability discussed in I Adams et al\ I I2007I) . 




Fig. 9. — The geometry of the toy model used to determine the torque on a test mass. The details of the model are described fully in 
Appendix 1X1 



